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Abstract. We study the stochastic evolution of four species in cyclic competition 
in a well mixed environment. In systems composed of a finite number N of particles 
these simple interaction rules result in a rich variety of extinction scenarios, from 
single species domination to coexistence between non-interacting species. Using exact 
results and numerical simulations we discuss the temporal evolution of the system for 
different values of N, for different values of the reaction rates, as well as for different 
initial conditions. As expected, the stochastic evolution is found to closely follow 
the mean-field result for large N, with notable deviations appearing in proximity of 
extinction events. Different ways of characterizing and predicting extinction events are 
discussed. 
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1. Introduction 

The study of evolutionary game theory and population dynamics often employs 
statistical physics and non-linear dynamics in order to shed light on multispecies 
ecological systems [U El EJ H]. By searching for generic characteristics in simplified 
models, we can start to understand coexistence and species extinction for complex, 
real-world systems. Many recent studies [51 EJ [H [HI El [10l HH [121 HSl US H51 [16l [TTl 
[IEl[19l[IBEni22l[23l[2^ have revealed 

a rich behavior in the evolution of three species in cyclic competition. This model can 
be expanded to produce more complex extinction scenarios by simply adding a fourth 
species J5j El EZl [38l EHl SHI S21 S31 SHI [36]. Labeling the different species with 
A,B,C,D, we allow A to 'prey on' B, B to 'prey on' C, etc. Interestingly, the four 
species form non-interacting partner-pairs, much like in the game of Bridge. Thus, in a 
system with N individuals, there are 2(iV + 1) absorbing states, most of which consist 
of a surviving partner-pair: A-C or B-D. This is of course fundamentally different to 
the three-species case where every species interacts with every other species. 

In [44] we used mean-field theory (MFT) to study the four species situation in 
absence of spatial and temporal fluctuations. MFT trajectories in the four-species case 
are influenced by a collective variable Q, see below, which evolves exponentially with 
an exponent A. A variety of trajectories in configuration space are encountered, ranging 
from periodic, saddle-shaped trajectories to spirals and even trajectories straight as an 
arrow. Although MFT is expected to capture some of the stochastic behavior, especially 
for large number N of particles, MFT cannot answer many interesting questions, 
and especially those related to extinction processes. In order to systematically study 
various extinction scenarios, we focus in the following on stochastic evolution by solving 
the master equation for small number of particles and by performing Monte Carlo 
simulations for larger systems. Varying predation rates, initial conditions, and the 
number of particles in the system allows us to develop a more comprehensive picture of 
the stochastic effects that take place in systems dominated by the formation of alliances 
of mutually neutral species. 

The paper is organized in the following way. We describe the details of the model in 
Section El As the update scheme for the stochastic processes is not unique, we propose 
two different schemes and elucidate the relationship between them. In Section [31 we 
summarize the mean-field theory results obtained in [H] as far as they are relevant 
for the present study. Section H] is the main part of the paper. We here describe in 
some detail different extinction events taking place in our system. Solving the master 
equation, we find exact results for a small system composed of N = 4 particles. For 
larger system sizes, we rely on Monte Carlo simulations to study the stochastic evolution. 
We study the behavior of Q during the process and discuss how systems may wind 
up on the 'wrong' — or unexpected — absorbing state. We also discuss two different 
maxims that are expected to predict the probable outcomes of our stochastic processes. 
We summarize our findings in the last Section and sketch possible avenues for future 
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research. 



2. The model and its absorbing states 

Our model involves four interacting species with cyclic competition. Denoting particles 
of each species by A,B,C, and D, the dynamics may consist of picking a random 
pair (PRP) and, if the individuals are " cyclically" different, letting them interact with 
probabilities p a , p^, etc.: 

A + B V A A + A 
B+C 4 B + B 
C+D^C+C 
D + A^ D + D 

We emphasize that AC and BD pairs are non-interacting. A mnemonic for these rules 
could be U A consumes B with probability p a \ B consumes C ..." No spatial structure 
is imposed, so the configuration of the system is given by the set of four integers 
(Na, N b , Nc, Nd) where Nx indicates the number of particles of type X present in 
the system. Note that the total number 

N = N A + N B + N C + N D (1) 

is an invariant, so that the configuration space is just a set of points within a three- 
dimensional simplex, namely a tetrahedron (of length N on each side). 

It is easy to write down the master equation for P (Na, N b , Nc, Nd',t), the 
probability for finding the system r steps after some given initial distribution, 

P (N A ,N B ,N C ,N D ): 

P(N A ,N B ,N c ,N D ;r + l) 
Pa (N A -1)(N B + 1) 



N (N — 1) /2 



■P(N a -1,N b + 1,N c ,N d ;t) 



+ MN B -l)(Nc + l) 
N(N-l)/2 
Pc (N c -l)(N D + l) DfAT 

N (N — 1) /2 ^ ' ' C ~ 1 > Nd + X ' t > 

+ N(N _ 1)/2 P(Na + 1,N b ,N c ,N d -1;t) 



+ 



Z 



N(N-l)/2 



P(N a ,N b ,N c ,N d ;t) (2) 



where 



Z = p a N A N B + p b N B N c + p c N c N D + p d N D N A ■ (3) 



From here, we can derive a partial differential equation for the generating function. But, 
to find its solution is far from trivial. As summarized in Section [3j applying a mean-field 
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approach has proven to be fruitful in providing insight into the behavior of the system, 

see nam]. 

A different approach is to perform computer simulations. However, it is clear that 
using the PRP scheme, many attempts to change the system would fail and the evolution 
would be quite slow. To speed up the process, we exploit another scheme, in which we 
always change a pair (ACP) at each step. Here, "time" t is measured in terms of 
interaction steps. To be precise, if the system at time t consists of Na (t), Nb (t), etc., 
we first construct Z(t), the combination above. We then generate a random number 
and, if it lies in the range 

p a N A (t)N B (t) \ 

Z(t) ) ' 1 ] 

set N a (t + 1) = Na (t) + 1 and Nb (t + 1) = Nb (t) — 1. Similarly, we increase and 
decrease, respectively, Nb and Nc, if it lies in the range 

' PoNa (t) N B (t) p a N A (t)N B (t) p b N B (t)N c (t) \ 

Z(t) ' Z(t) Z(t) J ' 1 j 

etc. In this manner, a configurational change occurs at each increment of t, with the 
appropriate relative probabilities. The master equation for P (Na, Nb, Nc, No',t), can 
be written for the ACP scheme in the following form: 

P(N A ,N B ,N c ,N D ;t + l) 
p a (N A - 1) (N B + 1 



Z+ Pa (N A -N B -l 



-P (N A -l,N B + l,Nc,N D -t) 



Z + p b (Nb -Nc-i) 

Z + p c (N c -N D -l) 

Pd (N D - 1) (N A + 1) , 

+ rk~r at ttP(N a + 1,Nb,N c ,N d -l;t) . 6 

Z + Pd (N D — Na — 1) 

However, the presence of variables in the denominators poses serious challenges, 
even for deriving an equation for the generating function. Needless to say, the details 
of the evolution under these two schemes will be very different. 

Let us remark that each face of the tetrahedron is "absorbing" in the sense that 
transitions into the face are irreversible and corresponds to the extinction of one 
of the four species. Within each face is a special limit of the problem with cyclic 
competition of three species, where one of the three rates is zero. It follows that lines 
between nodes indicate extinction of two species and vertices indicate the extinction 
of all but one species. Further, unlike the cyclic competition of three species, which 
only has three absorbing states regardless of N, our system has 2 (N + 1) absorbing 
states, namely, the two fixed lines: A-C and B-D. Therefore, it is a priori unclear if 
the transition probabilities, associated with an initial configuration ending up in any 
particular absorbing state, are the same for the two schemes. Fortunately, we are able to 
prove that these extinction probabilities are scheme-independent (see the Appendix for 
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details). Thus, we can use the fast, ACP scheme with confidence, although we should 
refrain from comparing every aspect of its full evolution with the time dependence from, 
say, a mean-field approximation of the PRP updates. 



3. Summary and implications of mean-field theory 

While writing the master equation (J2J) is facile, solving it or obtaining an appropriate 
generating function is very challenging. However, much insight into this problem can 
be gained by assuming a far simpler approach: the mean-field approximation, which 
focuses on the mean values of the fractions 

(X i ) r /iV = ^(Z i /iV)P(AS,C', J D;r) . (7) 

Mean field theory of our four species model was the topic of two recent publications. 
Therefore, we here only summarize the main results that are needed for the following 
discussion and refer the reader to [SI Hi] for more details. 

We find it convenient to simplify the notation by denoting these fractions by A(t), 
B(t), etc. from here on and alert the reader when necessity requires us to revert to 
integer values. Thus, we have 

A (r) + B (r) + C (r) + D (r) = 1 . (8) 

Following standard routes, we can derive an equation for the changes, (Xj) r+1 — 
(Xi) T , from equation ([2]). These will involve averages of products (e.g., (AB) T ) on the 
right. The mean-field approximation consists of neglecting all correlations and replacing 
the averages of products by the products of averages, so that the end result is a closed 
set of determinstic equations for the averages A(t), B(t), etc. Finally, by rescaling 
time with N and considering the N — > oo limit, r can be regarded as continuous and 
differences can be replaced by d T . The result is the mean- field approximation: 

d T A = [k a B - k d D] A (9) 

d T B = [k b C - k a A] B (10) 

d T C = [k c D - k b B] C (11) 

d T D = [k d A - k c C] D (12) 

where we suppressed the (r) in A(t), etc. In this form, the probabilities (fc's) can be 
thought of as rates, which are often normalized to k a + k b + k c + k d = 1 in the literature. 
These rates are related to the original probabilities via k m = 2p m [13]. 

Since exponential growth and decay in populations are common, it is natural to 
write equations (l9])- (jT2]) as 

d T In A = k a B - k d D; d T In C = k c D - k b B (13) 

d T In B = k b C - k a A; d T In D = k d A — k c C (14) 

which also exposes a coupling between the pairs AC and BD, which will often be 
referred to as partner pairs. Constructing appropriate linear combinations, we exploit 
an important control parameter, 

A = k a k c - k b k d (15) 
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which allows us to highlight the role that a single species has on the growth/decay of 
the opposing partner pair, producing 

d T [k b In A + k a In C] = XD; d T [k c \nA + k d In C] = XB 
d T [k c lnB + k b \nD] = - XA; d T [k d \nB + k a In D] = -AC. 

Adding and subtracting appropriately, we see that the quantity 

Q = Bk c +k d £)k a +k b 

evolves in an extremely simple manner: 

Q (r) = Q (0) e XT . (17) 

We believe Q is a generalization of the quantity R = A kb B kc C ka that has been 
introduced in [17] for the three species case, except that R is invariant for any set of 
rates. This is not the case with Q as it has a time dependence. Clearly, special properties 
will be manifested in the class with A = 0, and we will refer to them as "quasi-stationary 
systems." 

For those readers interested in deterministic trajectories, their origin, and in-depth 
analysis we point them to [43J. Since these results have been well- reported, we choose 
to only focus on the conclusions here. 



3.1. Periodic systems: A = 

For such systems, k a k c = kf,k d and consequently, the numerator and denominator of Q 
are both invariant. Thus, Q is an invariant as well. In the evolution of this system, 
MFT provides that none of the species go extinct and the system evolves periodically. 
Much like the three species case, there is a closed loop in phase space. The closed loop 
can be described by defining the constants of motion 

/ = A kb C ka ; g = B kd D ka , (18) 

which can be used to define hyperbolic sheets through the tetrahedron, and their 
intersection is the closed loop and resembles (the rim of) a saddle. 

To continue our characterization of such an orbit, we consider the extremal points 
of the time evolution of A, denoted by A±. At such a point, the values of the other 
three fractions are, in general, not extremal themselves. As the details have previously 
been reported we remark that the equation for fixing them is 

A kb/ka A ± + C Al kb/ka = (3 (19) 

where /3 is a constant and depends only on the rates and the initial conditions Bo and Dq. 
Once A± is known it is a simple procedure to obtain the remaining values. As A± varies 
from to 1 there are typically two solutions. When these solutions coincide, A + = A_, 
we are at the fixed line formed by the intersection of the two planes: k a A = k^C and 
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k a B = k^D. Exploiting A + B + C + D = 1, the fixed line is easily represented with the 
parameter 7 e [0, 1]: 

(^) = r]Tr = ^r(i-7). (20) 

Clearly, the line runs between the absorbing A — C and B — D lines. 

If the system begins in the neighborhood of the fixed line, the trajectory is nearly 
circular with angular frequency J k^k^ (1 — 7) and further away begins to represent the 
surface of the tetrahedron [36]. This provides the saddle shape as described earlier in 
typical systems, systems not in the neighborhood of the fixed line or near an absorbing 
boundary. 

We caution the reader that in this regime (A = 0) MFT is a good indicator of the 
behavior of a fully stochastic system only until an extinction event is near. The aim of 
Section H] is to address this issue. 



3.2. Spirals and arrows: A 7^ 

With such rates, Q grows/decays exponentially, so that there are no non-trivial fixed 
points. Since the fractions are bounded by unity, this implies that either BD or AC 
vanishes in the large t limit. In other words, for a system with finite number of 
individuals, extinction of one of the species must occur quite rapidly, see Figure 4a of 
[13] for an example. For many questions of interest, the mean-field approach faces, not 
surprisingly, serious limitations. For example, it cannot predict, given a finite system 
with stochastic dynamics, the (average) time for one of the four species to become 
extinct, or the time for the system to reach an absorbing state. Another interesting 
question is: At the time when one population goes extinct, what is the composition 
of the other three species. In other words, when the system "lands" on an absorbing 
face (of the tetrahedron), where does it land and what is the associated probability. 
Nevertheless, we can pursue its consequences once we are given where we land. For 
example, suppose D vanishes first and we land at (Ai,Bi,Ci) on the ABC face. Such 
a three-species system is much easier to analyze than those in [17] . since C does not 
consume A. In the mean-field approximation, the solution is trivial: From equation 
( TTBl . we see that A kb C ka is an invariant while A is monotonically increasing (C being 
the non-consumer). Since the evolution ends when B becomes extinct, the endpoint, 
(Af, 0, Cf) on the AC line is given by an equation for, say, Af : 

A h/k a Q_ Af > j = A kb/ka d . (21) 

This equation can be solved in general, graphically or numerically. Of the two solutions, 
the larger must be Af (since A is monotonically increasing). As will be shown below, 
these predictions are useful for a finite, stochastic system only when A is not too close 
to zero. 
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4. Extinction events 

The full stochastic problem of extinction probabilities can be stated as follows: Given 
a set of rates and an initial condition, what is the probability that the system will 
be found in a stationary final state? Furthermore, at what time is a species most 
likely to go extinct? And what is the average time for extinction? Needless to say, 
finding the exact answers to these questions is far from easy, even for the three-species 
model. In our case with four species, the problem is considerably more difficult, since 
there is a "macroscopically large" number (2N + 2) of absorbing states. For typical 
values of N in the hundreds or thousands, even exploring this issue with simulations 
is non-trivial, since the parameter space is six-dimensional for each TV (three for the 
rates and three for the initial fractions). It is therefore unrealistic to study the full 
parameter space. In this Section, we will provide insight into the stochastic behavior of a 
number of interesting cases, using both exact solutions and Monte Carlo simulations with 
ACP rates. Earlier publications jl3] and jS] reported two general maxims concerning 
extinction of a species: "The prey of the prey of the weakest is the least likely to survive" 
and "The prey of the prey of the strongest is quite likely to survive." In the following 
subsections, we further validate these claims, as well as describe their limitations. 

4-1. Exact results for N = 4 systems 

Let us first provide exact results for the simplest non-trivial system: N = 4, starting 
with one individual in each species: = Nb = Nc = Np = 1. We will not consider 
any other initial conditions, all of which contain only three species, with relatively trivial 
extinction probabilities. 

In this case, it is more convenient to use the populations, rather than fractions, 
to denote the system's configurations. Thus, the initial state is given by (1, 1, 1, 1). 
To simplify notation further, we set pd to unity and restrict the others rates to the 
unit cube: p a ,Pb,Pc £ (0,1]. This can be done without loss of generality, since cyclic 
symmetry allows us to label D as the biggest consumer. 

With N = 4, there are ten extinction probabilities. To compute these transition 
probabilities, we simply enumerate all trajectories and their associated weights. For 
example, the weight of (1, 0, 2, 1) ->■ (1, 0, 3, 0) is just 2p c / (2p c + p d ) = 2pJ (1 + 2p c ). 
Tabulated below are the exact transition probabilities from (1, 1, 1, 1) to the final states: 



(4,0,0,0) 



2Pc (3p a + Pb) Pb 



(0,4,0,0) 



cr (p a + 2p b ) (p a + p b ) (2p a + p b ) 
2(3p b + p c )p 2 c 



(0,0,4,0) 



a (p b + 2p c ) {j>b + p c ) (2p b + p c ) 
2p a (3 Pc + 1) 



(0,0,0,4) 



a (p c + 2) (p c + 1) (2p c + 1) 

2p 6 (3 +Pa)pl 



o (1 + 2p a ) (1 + p a ) (2+p a ) 
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(3,0,1,0) 
(2,0,2,0) PaP , 
(1,0,3,0) 
(0,3,0,1) 
(0,2,0,2) Pb 
(0,1,0,3) 



cr (p a + 2p b ) (p a + p b ) (2p a + p b ) 
p a + 2p b + p c + 2 



a ( Pc + 2) (p a + 2p b ) 



a (p c + 2) (p c + 1) (2p c + 1) 

k>\Vc 

cr (p b + 2p c ) (p b + p c ) (2p b + p c ) 
p b + 2p c + 1 + 2p a 
a(l + 2p a ) (p b + 2p c ) 



a(l + 2p a )(l+p a ) (2+p a ) 



where 



a = p a + p b + p c + 1 . (22) 

As an illustration, we consider the most symmetric case: all rates being unity. Then the 
transition probabilities to each of the vertices is 1/9; to the mid-points of the 2 fixed 
lines, 1/6; and to the rest of the points on the fixed lines, 1/18. 

Also instructive is the case with, say, p b <C 1 (i.e., B being the "weakest" and 
Pa,Pc = O (1)). Then, the final states and the associated probabilities are, to O (1), 



(0,4,0,0) 
(0,0,4,0) 



1 
a 



2p a (3 Pc + 1) 



(2,0,2,0) p c J - 



1,0,3,0) 



o-(p c + 2) (p c + l)(2p c + l) 

Pa + Pc + 2 



a (p c + 2) 



a (p c + 2) (p c + 1) (2p c + 1) 

A clear conclusion is that the "prey of its [the weakest, B] prey" has vanishing 
survival probability. Here, C is B's prey so that D is the "prey of B J s prey" - and goes 
extinct. Similar to the simple three-species case (where A survives with probability 
Pb/ (Pa +Pb+Pc)), we arrive at a intuitively reasonable maxim: The prey of the prey 
of the weakest is least likely to survive. By contrast, the 'law of stay-out' in [TT] 
("The species that is least frequently engaged in interactions has the highest chance 
to survive.") does not seem to always apply. In particular, if we let p a — Pb "C 1 (so that 
A is the least interactive species), the probability of B surviving is arguably higher. To 
be precise, the probability of all four survivors being B is 1/ (1 + p c ) to lowest order, 
whereas the probability of some survivors being A is p c j (1 +p c )- 
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C 



C 
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Figure 1. MFT evolution (thick blue line) and Monte Carlo simulation (dashed 
red line) for initial fractions (A,B,C,D) = (0.02,0.10,0.48,0.40) and rates 
(k a ,kb,k c ,k d ) = (0.4,0.4,0.1,0.1). Here, A = and so MFT traces out a closed, 
saddle-shaped orbit. The simulation initially follows MFT but then, due to stochastic 
noise, is driven to one of the absorbing faces, the system ending either with B — D 
coexistence (a) or with A — C coexistence (b). 



4-2. Stochastic simulations for A = 

Although it is futile to study the entire parameter space, even for A = 0, we can 
understand the general behavior of these unique systems by investigating a limited 
number of parameters. 

When A = 0, the mean-field approximation is quasi- stationary, yielding closed 
loops in the configuration tetrahedron. A stochastic trajectory will follow the closed 
loop at early times, as shown in Figure [U However, due to the intrinsic stochastic 
noise of the system, the trajectory is driven far from the mean-field trajectory. These 
quasi-periodic trajectories come to an abrupt stop when the trajectories hit one of the 
absorbing faces of the tetrahedron, resulting in the extinction of one species. In the 
resulting end game, a second species (the one not interacting with the extinct species), 
dies, yielding a final stationary state with two non-interacting species. During different 
trials of the same initial system, the noise happens randomly and changes the final 
distribution of species dramatically, as exemplified in Figure [1] where, despite starting 
with identical parameters, either (a) BD survives or (b) AC survives. Consequently, 
the idea of extinction probabilities does not make much sense with A = as extinction 
is based solely on the stochastic noise of the system and clearly the maxims found in 
earlier publications will not hold. 

4-3. Stochastic simulations for A ^ 

When A ^ 0, mean-field predicts that Q grows/decays exponentially Consequently, 
the stochastic system evolves quite rapidly towards an absorbing face and often ends in 
two-species coexistence. The extinction processes happen very early in the stochastic 
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Figure 2. (a) Mean field trajectory for the rates (k a , kb, k c , kj) — 
(0.1,0.0001,0.1,0.7999) and initial condition (A,B,C,D) = (0.1,0.7,0.1,0.1). The 
trajectory comes close to the ABC plane in the regions labeled 1, 2, and 3. (b) Time 
evolution of the fractions of the different species for the same case. The data are 
obtained using a fourth-order Runge-Kutta scheme with time step At = 10~ 5 . 



evolution, whereas in the mean-field approximation the system continues to evolve in 
a non-trivial fashion at later times. For this Section we choose to focus on extreme 
rates with an asymmetric initial condition. The following case is unique, in that the 
extinction processes only take place in very specific windows. This will allow us to 
highlight specifically the role that N — > oo has on the system, as well as to discuss 
interesting extinction scenarios. 

We consider a system with an initial population (A, B, C, D) = (0.1,0.7,0.1,0.1) 
and rates (k a , kb, k c , kd) = (0.1,0.0001,0.1,0.7999), resulting in a positive A. As 
expected, mean-field theory predicts the evolution of this system to spiral towards the 
A — C line, see Figure EK- In two distinct cases prior to the t — > oo limit, the trajectory 
comes very close to the ABC absorbing face (i.e. D's extinction). The time evolution 
of the fraction of the different species in mean-field approximation is shown in Figure 
|2b. In the time intervals when D is close to extinction, the fraction of species C is at 
a plateau. Once the D individuals recover, the number of C individuals also increases, 
which limits the number of Ds in the system and subsequently leads to another decrease 
of the number of Ds. 

It remains unclear to us how to use in the four-species case a methodology similar 
to that proposed in [12] for the three-species case that allows to determine how close 
the system comes to the absorbing face. As outlined in [12] we believe that for four 
species, the distance to the absorbing face is an extinction probability and scales with 
the number of competing individuals, N, but it is an open problem how to approach 
that problem analytically. 

Instead of an analytical approach, we therefore turn to Monte Carlo simulations in 
order to determine the behavior of the stochastic system. In these simulations extinction 
events exclusively take place in the time segments where the trajectory is near the ABC 
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0.1 0.2 0.3 0.1 0.2 0.3 0.1 0.2 0.3 
AAA 

Figure 3. Fractions of A and C particles the moment the D species becomes extinct, 
as obtained in Monte Carlo simulations of finite systems, with rates (k a ,kb,k c ,kd) = 
(0.1,0.0001,0.1,0.7999) and initial condition (A,B,C,D) = (0.1,0.7,0.1,0.1). The 
total number of particles are (a) N = 10 4 , (b) N = 10 5 , and (c) N = 10 6 . The clusters 
labeled 1, 2, and 3 correspond to the regions in Figure [5] that have the same labels. 
For each case 10000 independent runs were done. 

face, so that the fraction of D is very low and the fraction of C is constant. In Figure 
[3] we discuss three different cases characterized by different numbers of particles in the 
system: iV = 10 4 , 10 5 , 10 6 . For each system, we run 10000 independent simulations 
using the ACP scheme and plot the A — C composition at the moment when D becomes 
the first extinct species (the number of Bs follows directly from A + B + C — 1). In 
the smallest system (N = 10 4 ), we see three distinct clusters, denoted by 1, 2, and 3. 
The three clusters correspond to the three times the mean-field trajectory comes near 
the absorbing face, see Figure EJ Increasing the number of particles by 10, cluster 1 
completely vanishes and clusters 2 and 3 become better defined. Finally, in the N = 10 6 
case, only cluster 3 remains. It is this last cluster that evolves into the unique long-time 
mean-field solution in the limit N — > oo. The extinction events that lead to the two 
other clusters, however, are due exclusively to stochastic effects taking place in finite 
systems. As shown in Figure H] the probability P for a given run to end up in cluster 3 
varies as a function of system size N as P = 1 — e~ aN , with the numerically determined 
value a ~ 13.5 x 10 -6 . This result validates the claim that the distance to the absorbing 
face scales with N. 

Finally, let us briefly discuss how our system evolves once the D particles have 
died out. At this point, we are left with a three-species system where the number of 
A particles, which do not have anyone preying on them anymore, increases steadily. 
Concurrently, the population of B particles also decreases steadily, as they are eaten at 
a rather high rate by the As, but only reproduce with the very small rate In fact, 
this endgame is very well described by mean-field theory, and the time evolution of the 
A and C populations closely follows the mean-field scenario described in Section 13.21 
A kb C ka is approximately an invariant and the positions of the endpoints on the fixed 
line nicely agree with Equation ( 12 lj) . 
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Figure 4. Probability as a function of N that a system ends in the final cluster 
3. The rates are (k ai kb, k c , kd) = (0.1,0.0001,0.1,0.7999) and the initial condition is 
(A, B, C, D) = (0.1, 0.7, 0.1, 0.1). The line is a fit to the function P = 1 - e~ aN , where 
a 13.5 x 10~ 6 . 



4-4- The evolution of Q 

As already mentioned previously, the quantity 

J^k b +k c (jk d +k a 

Q = £k c +k d Dk a +k b ( 23 ) 

fully characterizes the fate of our system in mean-field approximation. Indeed, Q then 
grows or decays exponentially according to Q (r) = Q (0) e Al ~ where A is the combination 
of predation rates given in Equation ( 1151) . 

Before we can directly compare the stochastic evolution of Q to these mean-field 
results, we first need to scale the simulation time t to the mean-field time r. The 
rescaling is achieved by incrementing at a particular time t (the ACP time) the PRP 
time r by 

dr(t) = 1/Z(t) (24) 

such that 

T(t + l)=T(t) + dr(t). (25) 

Here Z(t) is given by the normalization ([3]). 

The difference between the evolution of Q(t) and Q{r) is striking, as shown in 
Figured! Only after rescaling time, do we find that ln(Q) evolves roughly linearly with 
a slope close to A. Of course, ln(Q) evolves more linearly for larger system sizes, as is 
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Figure 5. Evolution of ln(<5) versus (a) ACP time t and (b) PRP time r for a system 
composed of N = 10 5 particles, with rates (k a ,k b ,k c ,k d ) = (0.1,0.0001,0.1,0.7999) 
and initial condition (A, B,C, D) — (0.1,0.7,0.1,0.1). When plotted as a function of 
t, ln(Q) linearly increases until reaching positive infinity when D dies out (not shown). 
The final distribution of species for this run is 6% A and 94% C . 

expected from our earlier discussion on finite size effects. Despite this, Q itself is not 
a reliable indicator of full extinction scenarios. As soon as one species dies, Q becomes 
trivial but competition between the other three species continues. Other limitations of 
Q are exemplified in the following discussion. 

4-4- 1- Systems with A = 0. In mean-field theory, when A is zero, both the numerator 
and denominator of Q (and, therefore, Q itself) are invariant. However, in the stochastic 
process, when A is zero, Q wanders from its initial value due to stochastic noise. 
Figure [6] illustrates the spreading histogram of hi(Q(t)/Q(0)) as a function of the 
number of interactions. This distribution highlights the random fluctuations in ln(Q) 
that reflect the simulation trajectories wandering away from the closed saddle-shaped 
orbit predicted by mean-field theory. The spread in Figure M appears symmetrical, as 
expected for a purely random distribution. However, for a finite system this changes 
dramatically once extinction events take place. In fact, over 90% of the trials, that 
compose the distribution shown in Figure El ended with BD coexistence despite the 
seemingly unbiased spread of \a{Q)\ 

4-4- 2- Systems with A / 0. When the pair AC has a larger product of consumption 
rates than BD, A is positive and, consequently, Q grows exponentially, tending towards 
positive infinity as B or D approaches extinction. The opposite occurs when A is 
negative: Q exponentially decreases until reaching zero when A or C dies out. However, 
there are often exceptions to these generalities. For example, if A is positive, Q will 
exponentially increase; but, if A or C is the first species to die out, then Q will 
immediately go to zero and the simulation will abandon MFT. The same phenomenon 
occurs when A is negative and B or D dies first. Of course, the final value of Q depends 
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Figure 6. Probability distribution of \n(Q(t)/Q(0)) as a function of ACP time t. 1000 
independent runs were done, where A(0) — B(0) = C(0) = D(0) = 2500 and the rates 
are (0.4, 0.1, 0.1, 0.4), making A zero. 



on the population size N. Stochastic systems with small N are subject to finite size 
effects and are more likely to stray from MFT. 

A good example is provided by our "extreme rates" , with initial population fractions 
(0.1, 0.7, 0.1, 0.1). In that case, ln(Q) initially grows toward positive infinity, as expected 
for A > 0. However, for this initial condition it sometimes happens for small numbers 
of particles that the first extinction event is given by the dying of A, yielding B as 
the sole survivor. If that is the case, then ln(Q) jumps to negative infinity, its final 
value, despite its initial increase. In this sense, the stochastic evolution ended on the 
'wrong'-or unexpected-absorbing face. 



4-5. The maxims revisited 

As the parameter space for systems with three or more cyclically competing species 
is so huge, attempts have been made to summarize the possible outcomes by "laws" 
or "maxims." For example, for the three-species case a "law of the weakest" was 
proposed [17] that states that for the case of asymmetric interactions the "weakest" 
species (i.e. the predator with the smallest predation rate) survives with a probability 
that approaches one when the number of individuals gets large. This counter-intuitive 
observation is in fact due to the unique constellation of three species that all interact 
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with each other. Based on our earlier work [121 IB] , we proposed the following general 
maxim: "The prey of the prey of the weakest is the least likely to survive." Applying this 
maxim to the three-species case immediately yields the "law of the weakest." However, 
there are additional cases where the outcome of the simulations completely contradict 
the mean-field prediction. In order to deal with these cases, the following second maxim 
was proposed: "The prey of the prey of the strongest is most likely to survive." 

In order to get a better grasp at the validity of these maxims, we undertook a 
systematic study where we fixed the number of particles (N = 400) and the initial 
fractions (0.25,0.25,0.25,0.25), but varied the predation rates. We thereby always set 
the largest predation rate (that of species D) to kd — 1- In addition to the "extreme 
rates" discussed earlier, we considered the following values of A = k a k c — k b k d : —0.64, 
-0.16, -0.04, -0.01, 0, 0.01, 0.04, 0.16, and 0.64. For every value of A, multiple cases 
were studied, where we let the system run until a stationary state was reached. At least 
5000 independent trials were done for every case, and the survival probabilities were 
computed. 

Based on these data, we remark that the first maxim faithfully predicts the outcome 
for our small system when |A| is large, i.e. |A| = 0.16 or 0.64. The fate of the 
system in those cases is determined by A, and only in extremely rare occasions do 
we observe any deviations from the mean-field prediction. One such example is the case 
(k a , kb, k c , kd) = (0.1, 0.25, 0.9, 1) where in 3% of the trials B and D die out. From our 
previous discussion of finite-size effects, we expect that for larger system sizes less and 
less trials will end up in the "wrong" steady state. The second maxim also does very well 
for large negative A values, where it is by and large equivalent to the first maxim. For 
large positive A, however, the second maxim miserably fails. As the largest predation 
rate is k d = 1, large positive values of A mean that both k a and k c are rather large, 
whereas k b is very small. One example is given by (k a , k b , k c , k d ) = (0.7,0.19,0.5,1), 
yielding A = 0.16. As a result, A preys very efficiently on B whereas at the same time 
only few additional Bs result from the preying of B on C. Consequently, B has a very 
small probability to survive. 

For A close to zero, stochastic effects get more and more important, and there is an 
increasing probability that the system does not wind up in the stationary state predicted 
by mean-field theory. Consequently, the maxims get less reliable the closer to zero |A| is. 
Of course, for an increasing number of individuals we expect the first maxim to remain 
valid even close to |A| =0. 

5. Summary and Outlook 

We investigate the time dependent evolution and extinction probabilities of a simple 
model of population dynamics: N individuals of four different 'species', which compete 
cyclically. Though seemingly a trivial extension from a similar three-species game (often 
called rock-paper-scissors game), this system displays a much richer behavior. Since the 
four species form 'partner pairs', the stationary states consist of one of the pairs, with 
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N + 1 possible compositions in each case. As a result, there are 2(iV + 1) absorbing 
states with generally non-trivial distributions among them. In previous publications, 
we focused on mean-field theory to study the deterministic evolution of this nonlinear 
system. By manipulating the coupled differential equations, we found a single parameter 
A which controls the system's general behavior and the exponential evolution of Q. Of 
course, in MFT, no species ever goes extinct and therefore, we cannot study specific 
extinction scenarios. Thus, in this paper, we focused on stochastic models to explore 
extinction probabilities and timing. For small N, we solved the given Master Equation 
and gave the exact transition probabilities for one of these systems, where N = 4 and 
there is initially one individual of each species. 

For large N, we have to rely on Monte Carlo simulations. However, the full 
parameter space is seven-dimensional: three parameters for the predation rates (with 
k a + kb + k c + kd = 1), one for the population size N, and three for the initial population 
fractions (with A + B + C + D = 1). With so many variables, it is futile to explore 
the full parameter space, and we focused on a limited number of interesting cases. In 
many cases the systems evolve intuitively. In the limit iV — > oo, the stochastic evolution 
closely follows the mean-field prediction, as demonstrated in systems with "extreme 
rates" and N = 10 6 . Exploring in more detail the parameter space, we identified the 
following two maxims: "The prey of the prey of the weakest is the least likely to survive" 
and "The prey of the prey of the strongest is quite likely to survive." These maxims, 
however, do have limitations. For example, systems with A = follow the mean-field 
orbits at earlier times, but eventually wander away from the MFT trajectories and land 
on a random absorbing face due to the stochastic noise. These purely stochastic cases 
can not be described by maxims. 

Simply adding a fourth species to the popular three-species 'rock-paper-scissors' 
game reveals rich, complex nonlinear behavior and a tendency towards coexistence. The 
four-species case is characterized by the formation of two alliances composed of mutually 
neutral partners. Whereas we expect similar results for an even number of species 
that interact cyclically, thereby yielding two competing alliances, a more complicated 
situation prevails for an odd number of species. Among those cases, the three-species 
case is very special, as it is the only one where every species interacts with every other. 
It is therefore an interesting question what new phenomena emerge for five species, for 
example. Further interesting results can be obtained by putting these systems on one- 
or two-dimensional lattices. Beyond the straightforward situation of uniform interaction 
rates, one can imagine cases that are more representative of actual ecological systems 
where a prey actively moves away from its predators or where predators strategically 
chase their prey. Work along these directions is in preparation. 
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Appendix A. Equivalence of PRP and ACP schemes 

Let x = (A, B, C, D), with integer values, denote a configuration of our system. Starting 
with an initial configuration, x , we consider a trajectory in the ACP scheme that takes 
it to an absorbing state, x a , in n + 1 steps, through a sequence of configurations: 

T AC p : x , x^, x i2 , x in , x Q (A.l) 

By definition of ACP, the successive Xj's are distinct. The weight associated with this 
trajectory is 

X «l / V \ X «1 1 X «2 ) -P ( X in, ) (A-2) 

Here 

Pfex 3 ),^ (A.3) 

where k (x^ Xj) is the product of the rate (k) and the two appropriate populations 
associated with the change from Xj to Xj. Z (x«) is the sum of all possible transitions, 
i-e-, Ex, A;(x i ,x j ). So, for example, X; = (25,144,95,319) and x,- = (25,145,94,319), 
then k (Xi, Xj) = k b (144) (95) and Z ( Xi ) = k a (25) (144) + k b (144) (95) + k c{95) (319) + 
fcrf (319) (25). The overall transition probability, from x to x Q , is the sum of these 
weights over all such trajectories. 

Turning to PRP, we can associate a class of trajectories to each tacp 

r PRP '■ x 0i--- x 0) x ii)--- x iu x i 2 5 •••) X i»i-- x i„) x a (^-4) 

where x is repeated m times, x^ is repeated mi times, ... and x in is repeated m n 
times. In this class, each m ranges from to oo. The weight associated with such a 
trajectory is 

/ \m k ( x 0; x n) / \mii ( X »n X »2) / \m in ( x i ra ; x a) / a r \ 

^ So ^ iV(iV_i)/2 lSnj JV(JV-l)/2"" 1 <J N(N-l)/2 { ' 

where 

* s 1 - < A ^» 

is the probability for the system to stay at Xj. The next step is clear: summing over all 
trajectories in the class generates a factor of 

for each i. This factor combines with k (x,, x^) to produce p (x,, x^-), precisely the factor 
appearing in the associated t^cp- In other words, each trajectory in r PRP belongs to a 
class that can be associated with a unique tacp- In the limit of t — > oo, the sum (over 
this class) of their weights is exactly the weight for that tacp- Therefore, the transition 
probabilities (to go from any initial x to any particular final x a ) for PRP and ACP 
are identical. 
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